#include <iostream>
#include <cstdio>
#include <cstdlib>
//#include <cstdint>
#include <unistd.h>
#include <stdbool.h>
#include <string.h>
#include <math.h>
#include <omp.h>
#include <sys/time.h>
using namespace std;

#define PHI (24.0f/100.0f)
#define PI (3.14159f)
#define DEG2RAD(d,r) do{r=PI/180*d;}while(0)

#define KNRM "\x1B[0m"
#define KRED "\x1B[31m"
#define KGRN "\x1B[32m"
#define KYEL "\x1B[33m"
#define KBLU "\x1B[34m"
#define KMAG "\x1B[35m"
#define KCYN "\x1B[36m"
#define KWHT "\x1B[37m"

void usage(void)
{
	printf("Syntax Error!\n");
	printf("Usage: ./heatdist <iterations> <size> <heat radius> <heat value> <threads>\n");
	exit(-1);
}

int main(int argc, char *argv[])
{
	/* check programm call */
	if (6 != argc) usage();

	unsigned int i,j,k;
	double start, stop, runtime;
#if DEBUG == 1
	char color[50];
#endif
	
	unsigned int iterations = atoi(argv[1]);
	unsigned int viewsize = atoi(argv[2]);
	unsigned int heatspot = atoi(argv[3]);	// => m
	double heatval = atof(argv[4]);			// => H
	unsigned int threads = atoi(argv[5]);

	/* sanity check */
	if (iterations < 0 || iterations > 1e6) usage();
	if (viewsize < 0 || viewsize > 1e6) usage();
	if (heatspot < 0 || heatspot >= viewsize) usage();
	if (heatval < 0.0f || heatval > 127.0f) usage();
	if (threads < 1 || threads > 32) usage();

	/* adjust dimension to fit the desired size */
	unsigned int worksize = viewsize + 2;
	
	/* allocate memory */
	double **x, **xn, **swp, **circle;
	x = (double**) malloc(worksize * sizeof(double*));
	xn = (double**) malloc(worksize * sizeof(double*));
	circle = (double**) malloc(worksize * sizeof(double*));

	for (i=0; i<worksize; i++) 
		x[i] = (double*) malloc(worksize * sizeof(double));
	for (i=0; i<worksize; i++)
		xn[i] = (double*) malloc(worksize * sizeof(double));
	for (i=0; i<worksize; i++) 
		circle[i] = (double*) malloc(worksize * sizeof(double));

	/* init arrays */
	for (i=0; i<worksize; i++)
	{
		for (j=0; j<worksize; j++)
		{
			x[i][j] = 0.0f;
			xn[i][j] = 0.0f;
			circle[i][j] = 0.0f;
		}
	}

	/* set initial state of heat spot */
	int di, dj;
	double deg, rad;
	unsigned int origin;
	origin = worksize/2;
	for (i=0; i<heatspot; i++)
	{
		/* sweep all radii */
		for (deg = 0; deg <= 360; deg += 0.5f)
		{
			DEG2RAD(deg, rad);
			//printf("%lf %lf\n", deg, rad);
			di = origin + sin(rad)*i;
			dj = origin + cos(rad)*i;
			//printf("%i %i\n", di, dj);
			circle[di][dj] = heatval;
		}
	}

	omp_set_num_threads(threads);

	/* benchmark */
	start = omp_get_wtime();
	for (k=0; k<iterations; k++)
	{
		#pragma omp parallel private(i, j) shared (x, xn, circle)
		{
			unsigned int tid = omp_get_thread_num();
			unsigned int all = omp_get_num_threads();

			unsigned int start = (tid * worksize / all);
			unsigned int   end = ((tid+1) * worksize / all);

			if (tid == 0) start += 1;
			if (tid == all-1) end = worksize-1;

			/* keep heat source active */
			for (i=start; i<end; i++)
			{
				for (j=1; j<worksize-1; j++)
				{
					if (circle[i][j] != 0.0)
						x[i][j] = circle[i][j];
				}
			}
			#pragma omp barrier

			/* calculate next step */
			for (i=start; i<end; i++)
			{
				for (j=1; j<worksize-1; j++)
				{
					double *nv = &xn[i][j]; // next value
					double *cv = &x[i][j]; // current value
					double *cvd = &x[i+1][j]; // cv-down 
					double *cvu = &x[i-1][j]; // cv-up
					double *cvr = &x[i][j+1]; // cv-right
					double *cvl = &x[i][j-1]; // cv-left

					*nv = *cv + PHI * (-4.0f * *cv + *cvd + *cvu + *cvr + *cvl);
				}
			}
			#pragma omp barrier
		}
		
#if DEBUG == 1
		/* print output */
		for (i=1; i<worksize-1; i++)
		{
			for (j=1; j<worksize-1; j++)
			{
				double *cv = &x[i][j]; // current value

				if (*cv >= 0.0 && *cv < 18.0)
					strcpy(color, KWHT);
				else if (*cv > 18.0 && *cv < 36.0)
					strcpy(color, KCYN);
				else if (*cv > 36.0 && *cv < 54.0)
					strcpy(color, KMAG);
				else if (*cv > 54.0 && *cv < 72.0)
					strcpy(color, KBLU);
				else if (*cv > 72.0 && *cv < 90.0)
					strcpy(color, KGRN);
				else if (*cv > 90.0 && *cv < 108.0)
					strcpy(color, KYEL);
				else if (*cv > 108.0 && *cv <= 127.0)
					strcpy(color, KRED);
				else
					strcpy(color, KNRM);

				printf("%s%.3lf\t", color, x[i][j]);
				//printf("%s# ", color);
			}
			printf("\n\n");
			//printf("\n");
		}
		printf("iteration: %i/%i\n", k+1, iterations);
		usleep(1e3*50); // ms
		printf("\n\f");
#endif

		/* swap pointers for next iteration */
		swp = xn;
		xn = x;
		x = swp;
		swp = NULL;
	}
	stop = omp_get_wtime();

	/* free memory */
	for (i=0; i<worksize; i++)
	{
		free(x[i]);
		free(xn[i]);
		free(circle[i]);
	}
	free(x);
	free(xn);
	free(circle);

	runtime = stop - start;
	//printf("%i\t%i\t%lf\t%lf\t%i\n",
	//	viewsize, iterations, runtime, runtime/iterations, threads);

	std::cout	<< viewsize	<< "\t"
				<< threads	<< "\t"
				<< runtime	<< "\t"
				<< std::endl;
				

	return 0;
}

